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Abstract 

We apply recently developed smooth boundary conditions to the quantum 
Monte Carlo simulation of the two-dimensional Hubbard model. At half- 
filling, where there is no sign problem, we show that the thermodynamic limit 
is reached more rapidly with smooth rather than with periodic or open bound- 
ary conditions. Away from half-filling, where ordinarily the simulation cannot 
be carried out at low temperatures due to the existence of the sign problem, 
we show that smooth boundary conditions allow us to reach significantly lower 
temperatures. We examine pairing correlation functions away from half- filling 
in order to determine the possible existence of a superconducting state. On a 
10 X 10 lattice for [/ = 4, at a filling of (n) = 0.87 and an inverse temperature 
of /? = 10, we did find enhancement of the d-wave correlations with respect 

to the non-interacting possible sign of d-wave superconductivity. 

PACS Numbers: 05.30.Fk, 05.70.Fh, 71.10.-hx 
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I. INTRODUCTION 



The two-dimensional (2D) Hubbard model is considered to be one of the possible 
models to describe the high copper-oxide superconductors 0. Despite the fact that this 
model describes quite well the insulating state and some of the normal state properties of 
the copper-oxides, it is not clear whether it contains all the necessary microscopic features to 
also explain their superconducting properties 0. Although a variety of approximate calcu- 
lations 0-§| predict the existence of superconductivity in the doped Hubbard model, there 
is still need to find confirmation from a calculation which does not depend on uncontrolled 
approximations. 

Quantum Monte Carlo (QMC) simulations have emerged as a method of choice in trying 
to solve this model numerically, since, in principle, the interaction is treated in an exact 
way 0. The path integral formalism is the standard starting point, where the partition 
function is decomposed by using the Suzuki- Trotter formula [|]. A widely used algorithm 
employs the Hubbard-Stratonovich transformation to decouple the interaction and integrate 
out the fermionic fields , leading to a bosonic path integral with an effective action which 
depends on the fermion determinant given in terms of the bosonic fields only. The trace over 
the Stratonovich fields can then be replaced by a stochastic sampling in which the fermion 
determinant becomes essentially the weight of the distibution of Hubbard-Stratonovich fields 

However, these simulations have a number of hardware limitations that cannot be easily 
overcome even with the modern state-of-the-art supercomputers. The first problem is that 
they can only be implemented on relatively small lattices and at relatively high temperatures, 
since the computational time increases rapidly with the number of sites of the lattice, and 
as the temperature is lowered. Thus, ground state properties in the thermodynamic limit 
can only be obtained via finite-size scaling techniques and zero-temperature extrapolations. 
The second difficulty is the appearance of the sign problem when the system is doped away 
from half-filling |Tl[]. The sign problem arises from the fact that the fermion determinant. 
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which is interpreted as a probabhty weight of the distribution of Hubbard- Stratonovich 
fields, is not always positive definite, but decays exponentially with decreasing temperature. 
Thus, since averages must be computed by taking differences of two positive quantities of 
similar magnitude, the numerical problem becomes unstable, giving rise to large errors in the 
measurements. The sign problem in the Hubbard model becomes worse at low temperatures 
for values of the doping where a possible pairing phase is predicted to exist. 

The appearance of these problems has stimulated the search of new techniques for solving 
Hubbard-like models beyond quantum Monte Carlo simulations [0-0. In this paper we 
will instead show that some of the above limitations can be partially overcome using new 
types of boundary conditions (BCs) , smooth boundary conditions (SBCs) , which have been 
successfully applied to the study of one-dimensional (ID) systems within a large number of 
numerical techniques |15| . These BCs consist of smoothly decreasing the energy parameters 
appearing in the Hamiltonian as we approach the edges of the lattice. The result of this 
operation is that, instead of having a sharp and rigid boundary as is the case with open 
boundary conditions (OBCs), with SBCs the boundary extends itself into the system in such 
a way that its exact size is not fully determinable. In general, we will talk of the bulk of the 
system as the region where the energy parameters are constant, and of the boundary as the 
region over which the parameters are smoothly turned off. All measurements are made in the 
bulk region. In addition to reaching the thermodynamic limit on a relatively smaller system 
than with periodic boundary conditions (PBCs) or OBCs we find that with SBCs we obtain 
an improvement in the behavior of the sign problem which allows one to reach significantly 
lower temperatures and to explore the possible existence of a superconducting phase in the 
doped system. In Section II we introduce the Hubbard model with SBCs, in Section HI we 
demonstrate the validity of SBCs comparing them to PBCs in the discussion of the half-filled 
case, in Section IV we describe the behavior of the average sign of the simulation with SBCs 
and study pairing correlations in the nearly half-filled case, and, finally, in Section V we 
summarize and conclude. 
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II. THE MODEL WITH SMOOTH BOUNDARY CONDITIONS 



We consider the positive-f/ Hubbard model on a two-dimensional lattice defined by 
the Hamiltonian 

H = -Y. tiji4aCja + c]^Ci„) - ^/Xin^^ + ^ Ui{ni^ - ^){ni^ - (1) 

{ij),<J i 

which consists of a system of itinerant electrons with an on-site interaction of coupling 
strength f/j. Here is the nearest-neighbor hopping parameter between sites i and j. 
All the above energy parameters appearing in the Hamiltonian are scaled according to 
the smoothing function /j, shown in Fig. |1|, in such a way that Ui/U = fii/fi = fi and 
tij/t = \{,fi + /j), where t/, /i, and t are the bulk values. In the following, without loss 
of generality, we will take t = 1, and express all the other parameters appearing in the 
Hamiltonian of Eq.(|l]) in dimensionless form. The cl^. are fermion creation operators in a 
Wannier orbital centered at site i with spin cr, /i is the chemical potential, and njo- = clo-Qo- 
is the density operator. 

The above rescaling of the energy parameters in the Hamiltonian is motivated by a 



similar procedure that was successfully applied to several ID systems [T^. The choice of the 
smoothing function fi emerges from a generalization of the ID analogue, but is not unique. 
Here we use fi = y{l — di/[Nf + 1]), where di measures the distance of site i from the nearest 
edge of the system, starting at 1 for i on the outermost "frame" , or square, of sites. Nf is 
the number of frames in the boundary region. The smoothing function is given by |T^ 



y{x) 



1 X <= 



l + tanh§^] 0<a;<l (2) 



x>=l. 

The most important property of the the function y{x) is that it and all its derivatives are 
continuous everywhere. Thus in the limit of a large number of frames there are no precise 
interfaces or boundaries which would cause scattering. However, for most of the results 
shown here, we use Nf = 2, and fi takes on the values 0.817574476 and 0.182425524. A 



surprising result of our work is that one can obtain excellent smoothing using only two 
frames. 



III. THE HALF-FILLED CASE 

It is well known that the half-filled Hubbard model on an infinite lattice has an anti- 
ferromagnetic ground state with periodicity commensurate with the lattice ||16|. At finite 



temperature, from symmetry considerations of the antiferromagnetic order parameter, we 
know that the correlation length remains finite, and, thus, no phase transition occurs. On a 
finite lattice, however, the behavior is different, since the correlation length can only grow as 
large as the linear size of the system. Thus, we can still speak of an approximate transition 
temperature for a finite lattice corresponding to where the correlation length reaches the size 
of the lattice. This different behavior between a finite and infinite system is very crucial and 
could lead to the erroneous conclusion that even the infinite system has a finite transition 
temperature ||18|. In general, finite-size effects will be sensitive to the type of BCs used. 
Thus, it is important to choose BCs which elminate finite-size effects already on relatively 
small lattices so that the extrapolation to the infinite system can be taken without the need 
to consider extremely large lattices. 

Despite the fact that the behavior of the half-filled case is well known, we are interested in 
studying it in the context of SBCs since we want to show that we can obtain thermodynamic 
limit results on a relatively smaller lattice than when using PBCs. Also, since at half-filling 
there is no sign problem, we want to show that this effect is independent of the behavior of 
the average sign with SBCs relative to PBCs. We find that applying SBCs to a finite lattice 
indeed reduces finite-size effects and allows faster convergence to the thermodynamic limit. 
We consider two types of measurements, local quantities and correlation functions. We first 
consider the average kinetic energy per site, given by 

{Ki) = -^(E^*j(clc,. + (3) 
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where the sum runs over all j nearest neighbors of i. In Fig. |^ we show (Ki) as a function 
of the linear size of the system using PBCs and SBCs. The results for SBCs are obtained 
with the number of smoothing frames Nf fixed at 2. This means that the bulk will have 
a linear size that is 4 sites smaller than the corresponding system with PBCs. It is very 
striking that already on a 6 x 6 system with a 2 x 2 bulk region SBCs give a relatively good 
estimate to the value that we obtain on a 16 x 16 lattice with PBCs. A similar conclusion 
can be drawn from Fig. where we show the double occupancy {rii^nii). The finite-size 
effects in this case are much stronger with PBCs than SBCs. Notice that both in Figs. |^ and 
^ the measured quantities saturate to the same value for either type of BCs as the system 
size becomes larger, as expected, since the type of BCs should not have much effect on the 
behavior of a very large system. 

We also consider several other types of local measurements, leading to the same conclu- 
sion: the thermodynamic limit is reached on a smaller lattice with SBCs rather than with 
PBCs or OBCs. In Fig. ^ we show the average total energy as a function of the inverse 
temperature. Also in this case, we see that for each temperature considered, with SBCs we 
reach the thermodynamic limit more rapidly than with PBCs. This behavior is very similar 
to what we found in the numerical density-matrix renormalization group (DMRG) study 
of the spin-| Heisenberg chain with SBCs. Also in that case, with SBCs the ground-state 



energy is reached on smaller lattices than with OBCs ||17|| . 

We also examine various correlation functions in order to show that the behavior is 
similar to that of the local quantities. When using SBCs the correlation functions can only 
be evaluated for two sites separated by a distance only as large as the size of the bulk. 
Similarly, with PBCs we can only calculate correlations between points as far out as half the 
linear size of the system, since the lattice wraps around onto itself. An additional difference 
is that, since SBCs break the translational invariance of the Hamiltonian, the measurement 
of local quantities and correlation functions cannot be averaged over the entire lattice, as 
can be done with PBCs. In this case, the averaging can only be done over sites, or pairs 
of sites, which are equivalent under the reduced symmetry of the lattice (reflection and 



inversion) with SBCs. This imphes that, in order to obtain better statistics, we need to run 
the simulation longer than with PBCs. This is probably the only significant disadvantage 
of SBCs over PBCs. In Fig. |^ we show the spin-spin correlation function, defined by 

c{l) = c{l^, ly) = - n,+i^i){ni^^ - (4) 

for several systems sizes with PBCs and compare it to the case with SBCs on a 10 x 10 lattice 
with two smoothing frames. Again, it appears that with SBCs we can obtain on a smaller 
system the same results as with PBCs on a considerably larger system. Having shown that 
the use of SBCs allows one to reduce finite-size effects on relatively small lattices, we now 
study the effect of SBCs on the average sign by considering the Hubbard model away from 
half-filling. 

IV. THE DOPED SYSTEM 

As one dopes the system, the long-range antiferromagnetic order which is present in the 
ground state of the half-filled band is abruptly destroyed and short-range incommensurate 
magnetic correlations set in [|l^]. Here we are mainly interested in studying the possible 
existence of a superconducting phase away from half-filling. Recently, there has been an 
enormous effort in trying to obtain the phase diagram of the doped Hubbard model. Most 
of the studies have not been able to determine conclusively whether there is or there is not 
a superconducting phase. Approximate techniques have been used, such as self-consistent 
calculations, mean- field theories and conserving approximations ||12[, often suggesting the 
possibility of a parameter regime having d-wave superconductivity. QMC results, which 
show evidence of an attractive pairing interaction, but no direct evidence of pairing, have 
not been completely satisfactory because of the appearance of the sign problem already 
at relatively high temperatures. The possibility exists that the temperature at which the 
pairing occurs is significantly lower than the accessible temperatures even on as small lattices 
as an 8 X 8. 
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Although initially we did not expect changing boundary conditions would have a large 
effect on the average sign, in fact it can. In general, the average sign decays exponentially 
with the inverse temperature (3 and system size A^, {S) ~ e~^^^ |Tl||, where A depends on 
the QMC method (determinantal, world-line, projector, etc.). If {S) falls below about 0.1, it 
becomes impossible to obtain reasonable results. Previously, we showed that for large values 
of U the sign decays slower in the world-line method than in the determinantal approach 
20| . Here we find that the decay of the average sign can be substantially slower with SBCs 



than with either open BCs or with PBCs. 

In Fig. ^ we show the average sign as a function of inverse temperature for various lattice 
sizes at a filling of (n) = 0.87 with U = 4 for several types of BCs. If we consider the case 
with OBCs and PBCs only, we see that on an 8 x 8 lattice the average sign is already too 
small at a /3 ~ 6 — 7 for calculating measurements accurately. On the other hand, when we 
consider SBCs we see that the sign allows one to go to lower temperatures even on a lattice 
as large as a 12 x 12 with two smoothing frames. In general we find that the average sign as 
a function of temperature decays slower with SBCs than with PBCs and OBCs for all values 
of U. Particularly suprising is that the smooth system with a bulk region of size 4x4 (the 
solid triangles) has a sign that decays substantially more slowly than either a periodic or 
open system whose full size is 4 x 4. 

There are several ways in which one can understand why the average sign is better 
behaved with SBCs rather than with other types of BCs. When we decrease the energy 
parameters in the Hamiltonian on the edges of the system, we effectively make the local 
temperature there higher. Thus, smoothing the energy parameters of the system is equivalent 
to introducing a temperature gradient on the lattice without any heat flow into or out of 
the bulk of the system. (This strange temperature gradient cannot be realized in a real 
system, as far as we know.) Thus, if one thought of the average sign as a measurement 
that is obtained from averaging over the entire lattice, it might be reasonable to expect its 
behavior to be better with SBCs rather than PBCs. However, the average sign does not 
really involve some averaging process over the entire lattice — something slightly more subtle 

8 



is going on. 

However, this point of view can give us some insight into understanding why finite- 
size effects are much smaller with SBCs rather than with PBCs. Since the boundary is 
at a relatively higher temperature, a larger number of local states in the Hilbert space are 
accessible there. Thus, on the boundary, the system has the freedom to choose a larger 
portion of the Hilbert space, and thus, it can adjust itself in order to satisfy the constraints 
of the bulk. In the partition function, the lower temperature of the bulk region dominates 
the high temperature edge region, so the edges adjust to satisfy the bulk. 

As mentioned above, the most surprising effect that emerges from Fig. ^ is that the sign 
on an 8 X 8 system with SBCs and a 4 x 4 bulk region is better behaved than on a 4 x 4 
system with OBCs or PBCs. At this point we only have a qualitative argument for this 
effect. Based on the projector QMC simulation of the same system, we can at least show 
that this behavior is not completely unexpected. We have verified that the average sign in 
the projector QMC has a very similar behavior, namely that it is improved with SBCs. In 
the projector QMC the sign problem originates from the phase of matrix elements (entering 
in the partition function) of the form, 

(0(O)|0(r))^(0(O)|e--^|0(O)), (5) 

where \4>{0)) is an initial trial state (which can be, for example, a filled Fermi sea), and 
|0(r)) is the initial state propagated to imaginary time r = LAt. The fiuctuating Hubbard- 
Stratonovich fields cause the state |0(t)) to evolve through Hilbert space, and its precise 
evolution will depend on the type of BCs used. In particular, with SBCs, in the boundary 
region, where the temperature is high, the system evolves slowly, so that the propagated 
state is still very close to the initial state there. The state could evolve rapidly in the bulk 
region, except for the fact that it is continuously connected to the boundary region. Thus 
the edges act as a drag force on the bulk, slowing the rapid variations in imaginary time 
that cause the sign problem. 

In order to verify this picture, we studied, using the projector method, a simphfied "toy" 
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model: a 2D non-interacting tight-binding system with the addition of a Ising-hke field 
Aj = ±1 coupled to the density operator rii, designed to mimic the Hubbard- Stratonovich 



Two particles were put in the system, since a single particle never has a sign problem. 
The state of the higher energy particle had a nodal line which moved as the system evolved 
through imaginary time. If the nodal line rotated through 180°, the system had a minus sign. 
A specific field configuration was chosen with a spiral configuration through imaginary time 
to try to drag the particles in a circle in order to get a minus sign as quickly as possible. The 
field was able to cause rapid rotation of the nodal line in both PBCs and OBCs. In SBCs, 
however, the nodal line was unable to evolve rapidly on the edges (it was stuck, as if in a 
viscous fiuid), and the part of the nodal line in the bulk was held back. Consequently, it was 
much more difficult to generate configurations corresponding to minus signs. We believe the 
results support our picture for the improved behavior of the sign in the interacting system. 

We now present local measurements and correlation functions at high enough tempera- 
tures to be accessible to both PBCs and SBCs. When comparing the results with the two 
types of BCs we find similar behavior to the half-filled case. Then we will present results at 
lower temperatures, unaccessible with PBCs, and continue the analysis with SBCs. 

In Fig. we show the kinetic energy as a function of inverse temperature at a filling of 
(n) = 0.87, with U = 4. With PBCs we can reach only a temperature of /3 = 6 on an 8 x 8 
lattice, while with SBCs we can reach /3 = 12. 

We considered three types of pairing correlation functions, corresponding to different 
symmetries of the order parameter. In general, a given pairing correlation function is given 



field: 




(6) 



by 



D{1) = \{A,^iAl)\, 



(7) 



where the pairing field operators, A^, are given by 
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Af 



I - ci^cii, 



(8) 



for s-wave symmetry, and by 




(QtQ+^.| + Ci^Ci-^i) + - {ci^Ci+yi + Ci^Ci-yi) , 



(9) 



and 




(10) 



for extended s-wave and d-wave symmetry, respectively. 

In Fig. 1^ we show the pairing correlation functions for the above three different symmetry 
channels. For all temperatures considered, we find that pairing in a ci-wave channel is 
always an order of magnitude stronger than for an s-wave or extended s-wave channels. 
Thus, we continue our analysis with the rf-wave pairing correlation function only, examine 
its temperature dependence and compare it to the non interacting case to see whether 
there is an enhancement with respect to the U = case. In Fig. ^ we show the rf-wave 
pairing correlation function for several temperatures with U = 8, showing that, as the 
temperature is decreased there is an enhancement in the pairing. For reference, we also 
show the corresponding U = results to show that there is no enhancement relative to 
the non-interacting case, for the accessible temperatures. In Fig. we show the d-wave 
pairing correlation function for lower temperatures with U = A. Here, at (3 = 10, we do find 
evidence for enhancement over the non-interacting case, near the M point. This may be a 
sign of (i-wave superconductivity. 



We have implemented a numerical simulation of the 2D Hubbard model using SBCs. We 
have showed that at half-filling, where there is no sign problem, one obtains thermodynamic 
limit results on a smaller lattice than when using PBCs. Away from half-filling, we have 
found that the average sign decays more slowly with inverse temperature and lattice size with 



V. CONCLUSIONS 
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SBCs than with OBCs and PBCs, allowing us to reach significantly lower temperatures and 
larger lattices. We looked at the pairing correlation functions and showed that the d-wave 
channel is favored over other types of pairing channels, and that the pairing increases as we 
lower the temperature. On a 10 x 10 lattice for [/ = 4, at a inverse temperature oi P — 10, 
we did find enhancement of the d-wscve correlations with respect to the non-interacting case, 
a possible sign of d-wave superconductivity. 
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FIGURES 

FIG. 1. The smoothing function, /j = fi^^iy for a 12 x 12 lattice with four smoothing frames. 

The bulk region is a 4 x 4 lattice. The four concentric square frames are clearly visible. 

FIG. 2. The kinetic energy, (Ki), measured at the center of the lattice, as a function of linear 
system size with U = 4 at half-filling at /3 = 6. The solid squares are with PBCs and the solid 
circles with SBCs with two smoothing frames. 

FIG. 3. The double occupancy (nj-fn^), measuread at the center of the lattice, as a function 
of linear system size Nx with [/ = 4 at half-filling at /3 = 6. The solid squares are with PBCs and 
the solid circles with SBCs with two smoothing frames. 

FIG. 4. The total energy (-Ej), measuread at a site i in the center of the lattice, as a function 
of inverse temperature /3 for several system sizes and with PBCs and SBCs. The filling is (n) = 1 
and U = A. The open symbols are with PBCs and the solid symbols with SBCs. With SBCs we 
use two smoothing frames on the boundary. The solid lines are just guides to the eye. 

FIG. 5. The spin-spin correlation function c{lx,ly) at half-filling with U = 4 and at /3 = 6. 
The path taken in calculating c{lx,ly) is shown in the inset. The solid squares are for a 10 x 10 
system with SBCs and two smoothing frames. The open symbols are with PBCs and on sizes as 
shown. The solid line is a guide to the eye for the case with SBCs. The errors (not shown) are 
smaller than the size of the symbols. 

FIG. 6. The average sign {S) as a function of (3 with U = 4 and at (n) = 0.87 for different 
system sizes and BCs. The open symbols are for OBCs and PBCs and the solid symbols are for 
SBCs for the sizes indicated. The solid lines are just guides to the eye. 



15 



FIG. 7. The kinetic energy (i^i), measuread at the center of the lattice, as a function of /? 
for several systems sizes and with PBCs and SBCs. The filling is (n) = 0.87 and U = 4. The 
open symbols are with PBCs and the solid symbols with SBCs. The triangles are for an 8 x 8 and 
the upside-down triangles for a 10 x 10 lattices. With SBCs we use two smoothing frames on the 
boundary. 

FIG. 8. The pairing correlation functions D{1) with U = 4 at (n) = 0.87 and at /? = 8. The 
solid circles are for s-wave, the solid triangles are for extended s-wave, and the solid squares are 
for d-wave. The path taken through the lattice corresponds to a triangle as indicated in the inset. 
The corresponding empty symbols are for U = 0. The solid lines are just guides to the eye. 

FIG. 9. The d-wave pairing correlation function D'''{1) on a 10 x 10 lattice with U = 8 at 
(n) = 0.87 for several values of /3. The path taken through the lattice corresponds to a triangle as 
indicated in the inset. The corresponding empty symbols are for U = 0. 

FIG. 10. The d-wave pairing correlation function D'^{1) on a 10 x 10 lattice with ?7 = 4 at 
(n) = 0.87 for several values of The path taken through the lattice corresponds to a triangle as 
indicated in the inset. The solid lines are just guides to the eye. The solid squares are for /? = 4, 
the solid circles are for /3 = 6, the solid triangles are for /3 = 8 and the open squares are for P = 10. 
The corresponding empty symbols are for U = 0. 
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